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ABSTRACT 

We tighten previous upper limits on gamma-ray burst repetition by analyzing the angular power spec- 
trum of the BATSE 3B catalog of 1122 bursts. At 95% confidence, we find that no more than 2% of all 
observed bursts can be labeled as repeaters, even if no sources are observed to repeat more than once. If 
a fraction / of all observed bursts can be labeled as repeaters that are observed to burst v times each, 
then all models with (v — l)/> 0.05 are ruled out at 99% confidence, as compared to the best previous 
99% limit (v — 1)/ > 0.27. At 95% confidence, our new limit is (v — 1 )/> 0,02. Thus, even a cluster of six 
events from a single source would have caused excess power above that present in the 3B catalog. We 
conclude that the current BATSE data are consistent with no repetition of classical gamma-ray bursts 
and that any repeater model is severely constrained by the near-perfect isotropy of their angular dis- 
tribution. 

Subject headings: gamma rays: bursts — methods: statistical 


1. INTRODUCTION 

The origin of cosmic y-ray bursts (GRBs) is not known, 
and their distance scale has not been established. The 
angular isotropy of GRBs provides an important clue 
which has generated a “great debate” about the question 
whether GRBs are of Galactic or cosmological origin (e.g., 
Briggs 1995; Fishman 1995; Fishman & Meegan 1995; 
Hartmann 1995; Lamb 1995; Paczynski 1995). While 
cosmological models usually invoke singular events, such as 
the merger of two compact objects (e.g., Narayan, Paczyn- 
ski, & Piran 1992; Meszaros & Rees 1993), Galactic 
models currently under consideration require multiple out- 
bursts from each source. Recurrence in the framework of 
cosmological models could occur due to lensing, but the 
frequency of such events should be very small (e.g., Nemi- 
roff et al. 1994). Detection of a significant fraction of repeat- 
ing GRBs would argue against a cosmological origin and 
favor a Galactic (halo) origin. It would definitely exclude 
cosmological models in which the source is destroyed. 

To satisfy the isotropy constraint, a Galactic halo must 
be very large in order to minimize the dipole due to the 
solar offset from the Galactic center. The current multipole 
limits (Briggs et al. 1996; Tegmark et al. 1996) require 
Galactocentric shells with typical radii ^200 kpc. On the 
other hand, halos that are too large will yield an excess of 
bursts toward M31, which is not observed (e.g, Hakkila et 
al. 1994, 1995; Briggs et al. 1996). Because of these twin 
constraints, most halo models invoke a limiting sampling 
distance of about 300 kpc for the BATSE bursts. 

Currently, the only surviving Galactic model invokes 
bursts that arc produced by high-velocity pulsars (HVPs) 
born in the vicinity of the disk streaming out into the halo 


with velocities ~ 10 3 km s~ 1 (Li & Dermer 1992; Duncan & 
Thompson 1992; Duncan, Li, & Thompson 1993; Woosley 
1993; Li, Duncan, & Thompson 1994; Colgate & Leonard 
1994, 1996; Li & Duncan 1996; Bulik & Lamb 1996; Pod- 
siadlowski, Rees, & Ruderman 1995; Woosley & Herant 
1996). The recent upward revision of radio pulsar velocities 
(Lyne & Lorimer 1994) provides some support and motiva- 
tion for such a scenario. 

Radio pulsars are born in the Galaxy at the rate of 
roughly one pulsar every 100 years (e.g., Narayan & 
Ostriker 1990). This is consistent with recent estimates of 
the Galactic supernova rate [(2.5 ± 0.5) x 10' 2 yr 1 ; 
Tammann, Loffier, & Schroder 1994] if a significant frac- 
tion of them leave a black hole instead of a neutron star A 
typical value for the fraction of pulsars with velocities > I0 3 
km s' 1 is ^10% (Lyne & Lorimer 1994; Frail, Goss, & 
Whiteoak 1994). However, not every HVP may become a 
GRB source if additional selection criteria, such as particu- 
lar magnetic field strengths, must be applied. The fraction of 
pulsars with sufficient energy, if rotational momentum is the 
energy source, is probably much smaller than unity 
(Hartmann & Narayan 1996). On the other hand, BATSE 
has been detecting GRBs at the rate of approximately one 
per day. After correcting for Earth blockage and temporal 
gaps, the inferred all-sky GRB rate at the BATSE sensitivity 
level is ^ JO 3 yr" L The comparison of the pulsar birthrate 
w'ith the high observed rate of bursts implies that burst 
sources must repeat in the HVP scenario (Hartmann & 
Narayan 1996; Lamb 1995; Podsiadlowski et al. 1995). 
While the actual detection of burst recurrence is perhaps 
still avoidable, the parameter range of realistic halo models 
clearly encourages the search for repeaters in the data. On 
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the other hand, the leading cosmological scenario of 
merging compact objects would be in serious trouble if even 
a small number of GRBs could convincingly be shown to 
originate from the same source. 

Do classical GRBs repeat, and if yes, can we determine 
their recurrence pattern? The search for recurrence of clas- 
sical GRBs (in this work, we exclude the exciting class of 
soft gamma repeaters, which were recently reviewed by 
Kouveliotou 1994) has a long tradition (Mazets et al. 1981 ; 
Schaefer & Cline 1985; Atteia et al. 1985, 1987). Though 
somewhat model dependent, the distribution of burst loca- 
tions studied by these authors suggests a lower limit to the 
GRB recurrence time of — 10 yr. Interest in this topic was 
revived by Quashnock & Lamb (1993), who found evidence 
for repetition in the BATSE IB data set (Fishman et al. 
1994) using the nearest neighbor (NN) statistic (e.g., Scott & 
Tout 1989). 

Recurrence, even in a single case, would be immediately 
obvious if we had perfect locations. The locations provided 
by BATSE, while numerous, are imperfect, and consequent- 
ly a statistical analysis is required to demonstrate or limit 
the presence of repeaters. The positions in the first BATSE 
catalog had a minimum uncertainty of ~4° due to system- 
atics (Fishman et al. 1994). The repetition analysis of 
Quashnock & Lamb has been controversial. For example, 
Narayan & Piran (1993, 1994) used an apparent excess of 
burst pairs with ~180° angular separation (“antipodal 
bursts”) to argue against the repeater hypothesis derived 
from the small-angle NN excess, but Quashnock & Lamb 
(1994) argued that real physical (Galactic) anisotropies in 
the catalog are responsible for the positive antipodal corre- 
lations. 

In contrast to BATSE, COMPTEL and EGRET can 
localize bursts to — V , but the event rate of these detectors is 
very small. However, the recent near coincidence of 
COMPTEL bursts GRB 930704 and GRB 940301 (Kippen 
et al. 1995a, b) suggests repetition, because in 3 years of 
operation such a coincidence had only a 3% chance prob- 
ability. The most recent compilation of 27 GRB positions 
observed with COMPTEL (Kippen et al. 1996) did not 
yield another pair of coincident bursts, implying a reduced 
significance of the first pair. 

While most studies focused on projected GRB positions, 
it is clear that there could also be a clustering effect in the 
time domain. This aspect was investigated first by Wang & 
Lingenfelter (1994, 1995a, b), who also found evidence for 
repetitions. 

Given the importance of these findings for burst models, 
confirmation in subsequent samples is essential. Angular 
correlation function studies as well as nearest neighbor 
methods applied to the 2B sample did not confirm the 
earlier claims and instead found the data to be consistent 
with no repetition (Blumenthal, Hartman, & Linder 1994; 
Meegan et al. 1995). The small-scale excess was reduced, 
and the antipodal excess also went away (Hartmann et al. 
1994). However, the 2B data suffered from the problem of 
large data gaps in time due to the failure of the tape record- 
ers aboard the Compton Gamma Ray Observatory (CGRO). 
The resulting lower exposure to bursts obviously reduces 
our ability to detect recurrent events. Taking these effects 
into account, Meegan et al. (1995) constrain the total frac- 
tion of repeaters among the 585 observed 2B bursts to 
/ < 20%. This limit is based on simple repeater pairs in all 
cases. If each repeating source produces more than two 


detectable events, the corresponding limit would be tighter. 
The data also do not show evidence for enhanced clustering 
in time/space (Brainerd et al. 1995), although different sta- 
tistical measures, such as the Mantel-Haenszel test, provide 
marginal evidence for joint temporal and angular clustering 
with a timescale of 4- 5 days (Petrosian & Efron 1995). 

A comparison of the nearest neighbor method and the 
angular correlation function method shows that stronger 
constraints are obtained from the latter (Meegan et al. 1995; 
Brainerd 1996). These findings contradict earlier statements 
that the NN statistic is the superior method for studies of 
small-scale clustering (Lamb, Quashnock, & Graziani 
1994). Many authors have introduced new statistical tech- 
niques; most have concluded that there is either weak or no 
evidence for burst repetition (Hartmann et al. 1995; Bennett 
& Rhie 1996; Efron & Petrosian 1995; Hurley et al. 1994). 

To make further progress in the analysis of angular dis- 
tribution data, it is necessary to remove the two most 
important obstacles in this analysis; data gaps and poor 
positions. The most recent set of BATSE data (the 3B 
catalog; Meegan et al. 1996) is free of data gaps other than 
those induced by the South Atlantic Anomaly passages, and 
a major effort to improve the positioning algorithm 
(Pendleton et al. 1995) has reduced the systematic error to 
1?6. This sample of bursts thus provides a solid base for 
tests of the repeater hypothesis. Analyses of the 3B data 
using the standard tools of NN statistic and angular corre- 
lation functions (Meegan et al. 1996) confirmed the conclu- 
sions derived from the 2B data: classical GRBs have not 
been observed to repeat. The same conclusion was reached 
in a study of 3B data using the matched pair statistic 
(Bennett & Rhie 1996b). While supporting the null hypothe- 
sis, these studies did not significantly improve the limits on 
the repeater fraction. It is the purpose of this paper to intro- 
duce a new statistical method that provides significantly 
more statistical power than the standard tools. We apply 
the method to the 3B catalog and derive improved limits on 
the repeater fraction. 

The remainder of this paper is organized as follows: in § 2 
we present our repeater statistic, and in § 3 we apply it to 
the 3B data set and discuss the results. 

2. METHOD 

As in Meegan et al. (1995), we find it convenient to work 
with a two-parameter family of repeater models. These 
models are specified by the parameters / and v, where /is the 
fraction of all observed bursts that can be labeled as repea- 
ters and v is the average number of observed events per 
source observed to repeat. Our comparison of these models 
with the data will proceed as follows: 

1. We select a method of reducing an entire data set into 
a single number R, a number which is sensitive to the type 
of burst clustering that repetition would produce. We chose 
the sign so that the more evidence there is for repetition in 
the data set, the larger R will be. 

2. For fixed values of / and v, we compute the probability 
distribution of R by making Monte Carlo simulations of 
mock BATSE 3B catalogs. 

3. We compute the observed value of R, denoted R obs , 
from the real BATSE 3B catalog. 

4. Combining the results of the two previous steps, we 
obtain the function p(/ v), defined as the probability that 
R < R ob * given that the true parameter values are / and v. 
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For instance, if p(f, v) = 5%, then we would expect to 
observe as low an R-value as we did, by chance, only 5% of 
the time, and conclude that the model {f, v) is ruled out 
with 95% confidence. 

5. By repeating this analysis for a grid of points in 
parameter space and making a contour plot of p(f y v), we 
obtain our final results, shown in Figure 2. 

Clearly, the success of this approach depends crucially on 
the choice of the statistic R. Any choice whatsoever will of 
course give statistically valid results as long as the mock 
catalogs and the real data are treated in the same way, but 
poor choices of R will not allow good discrimination 
between repeating and nonrepeating models. The best 
choice of R is clearly that which allows us to rule out as 
many incorrect models as possible, i.e., in statistics jargon, 
that which gives our statistical test the maximum power. 
Many repeater statistics have already been employed in the 
literature, as discussed in § 1. Currently, the standard tests 
use the nearest neighbor test or the two-point correlation 
function. These are both good choices, as they are sensitive 
to small-scale clustering of the type that repetition pro- 
duces. Here we will use another choice, which we will argue 
is even better. 

2.1. The Total Power Statistic 

In Tegmark et al. (1996, hereafter THBM96), a method 
was presented for computing the angular power spectrum 
C i of gamma-ray bursts in the presence of the position 
errors of BATSE. It was found that in terms of the power 
spectrum, burst repetition has a very simple signature: the 
power at all multipoles / is increased by the same amount. 
Therefore, a logical measure of burst repetition would be 
the sum (or, apart from an irrelevant multiplicative 
constant, the average) of the power in all multipoles, R = 
Yj°= o C/, i.e., the total power. However, the position errors 
make the estimates of high multipoles very noisy, and in 
THBM96 it was found that the shot noise error bars on the 
power spectrum estimates explode for /> 70. To be useful, 
R should be fairly insensitive to noise, so we clearly want to 
give less weight to the C t with large error bars, i.e., with 
large /-values. With this in mind, we propose the following 
repeater statistic: 



As was shown in THBM96, the minimum variance estimate 
of the spherical harmonic coefficent a lm when faced with 
location uncertainties is 

| Y im(rW f)dQ . (2) 

Here x is the smoothed hurst map , plotted in THBM96, 
which is simply a sky map of all bursts smeared out by their 
position uncertainties: 

xtr) = £ (3) 

i= i 

where in the approximation of a Gaussian beam function. 


B k (r) 


1 

2n<?i 


exp 



(4) 


where 0 is the angle between r and r t , the position of burst k. 
Up to an irrelevant additive constant (the shot noise h tm 


discussed in THBM96), 1 d lm | 2 is a measure of the power C, 
Just as in THBM96, 

N; ff = Y, ex P + 0] > (5 

k= 1 

where o k is the uncertainty in the position of burst k anc 
gives the effective number of bursts that are well enough 
localized to contribute information about C t . In THBM96 
the error bars on C { where found to scale as 1/Nf ff , so the 
weights (N^ f /N) 2 in our definition of R have the desired 
property of suppressing the influence of the noisy high 
part of the power spectrum, since N? f -► 0 as / -> oo. 

2.2. A Faster Way to Compute the Total Power 

Since the quantity R is a measure of the total fluctuation 
power in the burst distribution, which is a rather natural 
quantity, one may ask if there is a simpler way of computing 
it which circumvents time-consuming calculations of large 
numbers of spherical harmonics Y lm . Fortunately, the 
answer to this question turns out to be yes. Substituting 
equation (2) into equation (1) and using the spherical har- 
monic identity 

oo i r 

X Z l*iJ 2 = |x(r)| 2 </Q i6 

I = 0 m = - 1 J 

(the spherical version of Parseval’s theorem), we obtain the 
simple and useful result 

R = J x(r) 2 d£l . [7 

Equation (7) provides an intuitive way of understanding 
how the statistic R works. If all bursts are well separated, sc 
that their respective Gaussians hardly overlap, then R will 
simply be a sum of separate contributions from each burst 
and will be independent of the exact burst positions. If clus 
tering is present, however, then we obtain “constructive 
interference” where two Gaussians overlap, and since the 
integrand in equation (7) is squared, we find a larger contri 
bution than if they did not overlap. Also, since the Gauss 
ians B k are normalized as probability distributions (they 
integrate to unity), their peak values are larger for well 
localized bursts. This means that overlaps conribute more 
to the sum if they involve a well-localized burst (with a 
small <j k ). 

2.3. A Still Faster Way 

We found that. An factors aside, the R -statistic is simply 
the mean squared amplitude of the smoothed burst map 
Since we wish to make many thousands of Monte Carle 
simulations to obtain accurate estimates of the probability 
distributions of step 2, it is desirable to further accelerate 
the procedure of computing R. Substituting equation (3 
into equation (7), we obtain 

R = £ £ | BtfiBfiDdQ . <8 

»=i j=i J 

In the approximation that o k 4 1 rad % 60° (which is quite 
an accurate approximation, as typical values are a few 
degrees), the integral reduces to that of the product of twi 
Gaussians in the flat two-dimensional plane and can bt 
done analytically. This leaves us with the handy result 
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where 0 ti = cos' 1 (r, * f j), i.e., the angle between the two 
bursts. Finally, since replacing our statistic R by a mono- 
tonic function f(R) (for instance, rescaling R and subtrac- 
ting off a constant) will in no way change its ability to 
discriminate between models, let us redefine R to give 
numbers of a convenient magnitude. If there is no clustering 
and the total area of covered by all the error circles is 
so small that substantial overlap is unlikely, then the sum in 
equation (9) will be dominated by the “diagonal 11 terms 
with i = j and reduce to R % (l/47c) of 2 , an expression 
which is completely independent of the burst locations 
Dividing by this quantity and disposing of uninteresting 
additive and multiplicative constants, we thus redefine our 
R-statistic as 




IJL I/ 

2 K 2 + <rj)]l 


(of + <rj) . 


( 10 ) 


This is the expression that we use in our calculations. Apart 
from an irrelevant factor of 4, it is merely the ratio between 
the off-diagonal (i ^ j) and diagonal (i = j) elements in the 
sum in equation (9). It will clearly always be nonnegative, 
and in the absence of clustering, it will approach zero if the 
location errors o x do. 

The position uncertainties A 6 quoted in the BATSE 3B 
catalog are defined as the radius of the 1 o circle, i.e., of the 
circle that contains erf (1/^/2)^ 68% of the probability. 
Thus, in the limit o k < 1, the conversion between A0 and o 
is 


^={-2ln[l-erf(-L)J'%0,6. (11) 

Note that the values of A 6 quoted in the BATSE 3B catalog 
do not include the systematic error contribution of 1?6, 
which is to be added to the quoted values in quadrature. 

In defining our statistic R above, we have omitted a few 
elements that were used in THBM96, for instance the spa- 
tially varying exposure function h and the Fisher beam 
function. It should be emphasized that despite these omis- 
sions and the various approximations made (that o k <? 60' , 
etc.), our statistical statements will be 100% exact. This is 
because, as mentioned above, we are free to define the sta- 
tistic R however we want, as long as we make no approx- 
imations when generating the Monte Carlo catalogs and 
compute R in exactly the same way from these and from the 
real data. The only real constraint is that if we make R 
depart too much from the exact measure of the total power, 
it may no longer be as good a measure of clustering, and its 
ability to reject incorrect models will be weakened. 

2.4. The M ock C a (a logs 

According to the model, the N = 1122 bursts were 
caused by (1 — J')N nonrepeating objects (rounded to the 
nearest integer) and fN/v repeating objects that burst v 
times each. Therefore, we generate the mock BATSE 3B 
catalogs as follows: 

1. The (1 —f)N + J'N/v objects are distributed randomly 
across the sky, in a completely uncorrelated fashion, but 
with the point density modulated by the exposure function 
h{r) as described below. 

2. Each nonrepeating object is assigned one burst, and 
each repeating object is assigned v bursts. 


3. The 1122 position errors o k from the BATSE 3B data 
set are resorted in a random order, and one is assigned to 
each burst. 

4. Each burst is displaced from its true position by a 
random amount drawn from the probability distribution B k 
of equation (4). (In fact, we use the more accurate Fisher 
beam function described in THBM96, but it is virtually 
identical for o < 60°.) 

We truncate the number of bursts from the last repeating 
source so that each mock catalog contains exactly 1122 
bursts. 

2.5. The Exposure Function 

The sky exposure of BATSE is not quite uniform 
(Fishman et al. 1994). The BATSE experiment does not 
exclude any area of the sky, but due to blocking by the 
Earth and detector gaps during passages of the so-called 
South Atlantic Anomaly, some positions on the sky have a 
reduced probability for burst detection. We quantify this by 
the function n(r), the exposure function , defined as the 
expected number of bursts per steradian. It is simply pro- 
portional to the total exposure time that each patch of sky 
has received. It is known a priori and is not obtained from 
the observed burst distribution. 

Because of problems due to the loss of the spacecreaft 
tape recorders, the absolute efficiency has not been deter- 
mined since the release of the IB data set. However, the 
shape of the exposure function h is essentially independent 
of time, and since the shape is all that matters for the 
present analysis, we employ the IB estimate (Fishman et al. 
1994). This function h depends on declination only and is 
independent of right ascension. This means that in equato- 
rial coordinates, the multipole coefficients h lm vanish except 
when m = 0. The dominant deviation from uniformity is a 
quadrupole (n 2 o/^oo ~ 8.8%) depletion of bursts near the 
equator due to the shadowing of the sky by the Earth. The 
second largest anisotropy is a dipole moment (n 10 /fi 00 ^ 
4.5%) toward the Earth’s north pole, due to the South 
Atlantic Anomaly, which requires disabling triggers. Com- 
pared to the shot noise, the higher multipoles (/ > 3) are 
negligible {a l0 /a 00 < 1%), but for completeness, they have 
nonetheless been included in our analysis. 

We incorporate the effect of variable exposure into our 
mock surveys as follows: 

1. We compute the maximum value of the function fi(r) 
and denote it n max . 

2. Whenever a burst has been generated at a position r, 
we accept it with a probability P = h{r)/h max and reject it 
with a probability 1 - P. In other words, we generate a 
uniformly distributed random number u e [0, 1], and if u > 
h{r)/h max , then we generate a new burst position r and a new 
random number w, repeating the procedure until we obtain 
u < h(r)/« max . 

The net result is that, when averaging over many mock 
catalogs, there are on average h(r)dft bursts in a solid angle 
dQ. around r. 


3. RESULTS AND DISCUSSION 

Figure 1 shows the cumulative probability distribution of 
the /^-statistic for a range of models with v = 2. For each of 
these models, we generated 10 4 Monte Carlo catalogs from 
which we computed the R-statistic, producing a single curve 
in the figure. These cumulative probability distributions 
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Fig. 1. — Monte Carlo results. The cumulative probability distribution 
of our test statistic R is shown, as determined from Monte Carlo simula- 
tions, for a range of repeater models. The models have a fraction / of the 
bursts caused by repeating sources that are seen to burst twice (v = 2), 
where, from left to right, /= 0, 2%, 5%, 10%, 15%, 20%, and 25%. The 
hatched region contains the ^-values smaller than that observed in the 
BATSE 3B data set, so the probability that a model is consistent with the 
data can be read off as the intersection of its curve with the vertical line. We 
see that / = 2% is ruled out at approximately 95% confidence, since the 
intersection takes place near the horizontal 5% line. (Note that since these 
curves depend on N and <r k , new Monte Carlo simulations must be made 
to analyze a different data set.) 

F(RJ show the fraction of the R-values that are smaller 
than any one constant R*, so F(R J is simply the integral 
from zero to R* of the probability distribution for R. For 
instance, the median R-value is the point R + at which 
F(RJ = 0.5. 

The first thing to notice about Figure 1 is that R behaves 
as expected: as the repeater fraction /'increases, the distribu- 
tion of R-values shifts further to the right. Second, the value 
extracted from the real data, R obs ^ 0.3085, lies far to the 
left in the figure, which means that the data contain no 
evidence whatsoever for repetition. F(R obs ) % 5% for the 
model v = 2, / = 0.02, which means that this model pro- 
duces such a low R-value only 5% of the time, i.e., that this 
model is ruled out at 95% confidence. Models with higher 
/-values are of course even less consistent with the data. 

At first sight, it may appear disturbing that R obs is lower 
than the typical values obtained for the null model / = 0 
which has no repetition at all. Even if there are no repeaters, 
R will only be as low as 0.3085 a mere 12% of the time. 
Since no form of clustering could explain this (rather, a 
contrived model with some form of “anticlustering” would 
be needed), one is led to ask how unlikely it is that this 
happened merely by chance. The answer is, of course, that if 
the no-repetition hypothesis is true, then the probability of 
finding an R-value this far out in one of the tails of the 
distribution is 2 x 12%, i.e., it would happen about a 
quarter of the time and thus should not be a source of 
concern. A similar conclusion was reached by Hartmann & 
Epstein (1989), who also found an unusually isotropic burst 
sample in their analysis of Interplanetary Network (IPN) 


positions from Atteia et al. (1987). In contrast, since wi 
know a priori that repetition causes high rather than low 
R-values, we should use a single-sided (rather than double- 
sided) test when ruling out repeater models. 

The results of repeating the Monte Carlo procedure foi a 
two-dimensional grid of points in the (/ v) parameter space 
are summarized in Figure 2. In agreement with Meegan et 
al. (1995), the contours are seen to be fitted accurately by 
the hyperbolic profile (v — 1) oc 1// At 99% confidence, we 
have 

(v- l)/< 0.049 , (12, 

and at 95% confidence, we obtain (v — 1 )/< 0.018. Foi 
comparison, the contour of models ruled out at 99% con- 
fidence with the correlation function method applied to the 
same data set is shown, (v — 1 )/= 0.27 (Meegan et al. 1996) 
It is seen that the method presented here tightens signifi- 
cantly the constraints on any observable repeater popu- 
lation. While no statistical analysis can disprove the 
possibility of a very low repetition rate, the constraints on 
simple repeater models discussed here are now so severe 
that either recurrence is very rare among classical GRBs (in 
contrast to what the earlier studies suggested), or it has a 
very special pattern we can no longer find in the 3B data 
although we did in IB. The latter interpretation appears to 
be very contrived as long as we do not find theoretical 
support for such a special pattern. 

The straight lines in Figure 2 correspond to models with 
a fixed number (1, 2, 4, 8, and 16) of repeating sources in the 
data set. For instance, the leftmost line corresponds to the 
case in which all the repeated bursts are due to a single 
source, so that the repetition would manifest itself as a 
single cluster of v bursts in the BATSE catalog. The fact that 
this line intersects the 99% contour below v = 9 means that 
our constraints are now so strong that not even a single 
repeater with v = 9 is allowed. Similarly, even a single 
sixfold repeater is excluded at 95% confidence. 

3.1. How Robust Are the Results ? 

How robust are these results to changing various 
assumptions about the data? Playing the devil's advocate, 
are there any types of errors that could produce artificially 
strong constraints? Most types of data problems, for 
instance, neglected modulations in the exposure function, 
would have the opposite effect and give weaker upper limits 
on repetition, since they would create extra clustering in the 
data. There are basically only two exceptions: 

1. If the true exposure function h is more smooth than 
assumed; 

2. If the location errors are underestimated. 

To quantify the effects of the first possibility, we reran the 
analysis with h constant. This is, of course, the most extreme 
case possible, since there is nothing more uniform than a 
uniform distribution, and quite contrived, since we know 
that the South Atlantic Anomaly and Earth shadowing 
must have some modulating effects. Nonetheless, the 
changes were quite marginal, with the weakened limits 
being (v — 1 )/< 0.064 at 99% confidence. In other words, a 
problem of type 1 could at most weaken the limits from 
0.049 to 0.064, i.e., by about 30%. 

A more realistic source of concern is the second possi- 
bility. Graziani & Lamb (1995) analyzed the distribution of 
3B locations in comparison to the positions derived from 
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Fig. 2. — Excluded repeater models. The leftmost shaded region of parameter space shows the repeater models that are ruled out by the BATSE 3B data at 
95% confidence. The middle shaded region is excluded at 99% confidence. The shaded region to the upper right is that ruled out at 99% confidence using the 
correlation function (Meegan et aL 1996) instead of our jR-statistic. From left to right, the straight lines correspond to models with 1, 2, 4, 8, and 16 repeating 
sources in the data set, respectively. 


the IPN 3 network and conclude that the systematic error 
A0 O of the 3B data should be ~4° instead of the advertised 
1?6. In addition, there may be correlations in the data that 
suggest a brightness dependence for the systematic error, 
instead of the constant value suggested by Meegan et al. 
(1996). The studies by Graziani & Lamb (1996) did not take 
systematic effects in the IPN localization method into 
account and also do not incorporate the fact that some IPN 
locations are based upon the earlier BATSE 2B locations 
and thus may be biased against the 3B locations. Although 
we tend to agree more with the error budget prescribed by 
the BATSE Team, this is an unsettled question, and we have 
therefore studied the effect of larger errors on our repeater 
limits. For simplicity, just to estimate the magnitude of the 
effects, we maintained a constant systematic error A0 O (as 
suggested for the 3B data) and repeated the entire analysis 
for A0 O = 2°, 3°, 4°, and 5°. To be maximally conservative, 
we took n constant as well, and we found that for v = 2, the 
99% upper limit on the repeater fraction scaled approx- 
imately as 

/ < O.O6(A0 o /1.6)° 4 . (13) 

In other words, even if the true systematic location errors 
were as high as 5°, our limits would only be weakened by 
less than a factor of 2. 

3.2. Conclusions 

In summary, we have sharpened previous limits on GRB 
repetition by analyzing the improved BATSE positions of 


the 3B catalog (Meegan et al. 1996) with a new statistic 
based on the angular power spectrum. This method is more 
powerful than the NN statistic or correlation functions 
because it is in a sense a global method (e.g., Hartmann et 
al. 1995). The presence of clustering (on any scale) some- 
where on the sky reduces the density of sources everywhere 
on the sky (relative to an isotropic distribution). Any 
method that utilizes the impact of clustering on all angular 
scales is more sensitive to clustering than local methods, 
such as NN statistics that only measure neighbor excess 
very close to a given source, or angular correlations. The 
appearance of clusters other than those occurring by chance 
quickly leaves a mark on the overall, global angular power 
spectrum. Using this effect, we find an amazing isotropy of 
bursts, which is hard to satisfy with any model other than 
that of nonrepeating sources, providing strong evidence 
against the Galactic models currently under consideration. 
It was believed that spherical harmonic expansion would 
not provide a good tool for repeater studies (Lamb et al. 
1994) because the power is spread over many high harmo- 
nics. We emphasize that the R-statistic introduced here 
sums over power in all modes, and the Monte Carlo simula- 
tions demonstrate clearly that this approach in fact does 
provide a very powerful tool for clustering studies. 

As shown in Figure 2, the bulk of the previously allowed 
parameter space is now ruled out, and the new constraints 
are so tight that 95% confidence, no more than 2%/2 = 1% 
of the burst sources can have repeated in the data set, 
assuming that repeaters have a brightness distribution 
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typical of the entire catalog. Continued burst observations 
would allow further improvements in studies of this kind, 
but to really advance the field we must develop new mis- 
sions that provide much better GRB localizations. While 
the forthcoming High-Energy Transient Experiment mis- 
sion (e.g., Ricker et al. 1992) may provide very accurate 
localizations, the event rate of HETE will probably be too 
low for recurrence to be detected. NASA has initiated 
several concept studies for new GRB missions, most recog- 
nizing the need for accurate positions. Future GRB obser- 


vations may free us from the need to use statistical method 
to answer such basic questions as burst recurrence, but fo 
now we have no choice. The angular power method pre 
sented here should serve us well until the launch of a nev 
generation of GRB experiments. 
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